package fr.inserm.um915.bdibd.model;

import java.util.ArrayList;
import java.util.List;
import java.util.SortedMap;
import java.util.TreeMap;
import java.util.logging.Logger;

import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.distribution.BinomialDistributionImpl;
import org.apache.commons.math.distribution.ContinuousDistribution;
import org.apache.commons.math.distribution.Distribution;
import org.apache.commons.math.distribution.IntegerDistribution;
import org.apache.commons.math.distribution.NormalDistributionImpl;


/**
 * 
 * DataExtractorFactoryIKinship
 *
 */
public class DataExtractorFactoryIKinship
	extends DataExtractorFactory
	{
	@SuppressWarnings("unused")
	private static final Logger LOG=Logger.getLogger("fr.inserm.um915.bdibd");
	private boolean use_zScoreAA=true;
	
	
	 public static class Stat
	  	{
		private double zScoreAA=Double.NaN;
		private double zScore=Double.NaN;
		
		public boolean isNan()
			{
			
			return Double.isNaN(zScoreAA) || Double.isNaN(zScore);
			}
		
		public double lodScoreAA()
			{
			if(isNan()) return Double.NaN;
			return (this.zScoreAA*this.zScoreAA)/2.0*Math.log(10);
			}
		
		public double lodScore()
			{
			if(isNan()) return Double.NaN;
			return (this.zScore*this.zScore)/2.0*Math.log(10);
			}
	  	}
	
	private abstract class DataExtractorKinship
		extends DataExtractorFactory.AbstractDataExtractor
		{
		private SortedMap<Integer,IBDStatus> pairid2theoritical=new TreeMap<Integer, IBDStatus>();
		private TreeMap<Float,Integer> typesAA = new TreeMap<Float, Integer>();
		private TreeMap<Float,Integer> typesAU = new TreeMap<Float, Integer>();
		private SortedMap<Integer,PairOfIndividual> pairid2couple=new TreeMap<Integer, PairOfIndividual>();
		private int nbPairsAA=0;
		private int nbPairsAU=0;
		
	
		
		private int getIbd0(LinkageRow linkageRow,int pair_id)
			{
			return linkageRow.get(pair_id).getIbd0() < getTreshold() ? 0:1;
			}
		
		@Override
		public void setPairOfIndividualIds(List<Integer> pair_ids)
			{
			List<Integer> list_of_pair_ids=new ArrayList<Integer>();
			this.pairid2theoritical.clear();
			this.typesAA.clear();
			this.typesAU.clear();
			this.nbPairsAA=0;
			this.nbPairsAU=0;
			for(Integer pair_id:pair_ids)
				{
				IBDStatus status= getModel().getTheoriticalIBDStatusByPairId(getTransaction(), pair_id);
				if(status==null) throw new RuntimeException("No theo status for pair_id:"+pair_id);
				//discard those having a fixed probability
				if(status.getIbd0()==1f || status.getIbd1()==1f || status.getIbd2()==1f)
					{
					continue;
					}
				
				this.pairid2theoritical.put(pair_id,status);
				
				PairOfIndividual p= getModel().getPairOfIndividualByPairId(getTransaction(), pair_id);
				if(p==null) throw new RuntimeException("No pair of indis for pair_id:"+pair_id);
				this.pairid2couple.put(pair_id, p);
				list_of_pair_ids.add(pair_id);
				
				switch(p.getConcordance())
					{
					case AFFECTED:
						{
						Integer f= this.typesAA.get(status.getIbd0());
						this.typesAA.put(status.getIbd0(),1+(f==null?0:f));
						this.nbPairsAA++;
						break;
						}
					case DISCORDANT:
						{
						Integer f= this.typesAU.get(status.getIbd0());
						this.typesAU.put(status.getIbd0(),1+(f==null?0:f));
						this.nbPairsAU++;
						break;
						}
					}
				
				super.setPairOfIndividualIds(list_of_pair_ids);
				}
			}
		
		
		public Stat getStats(LinkageRow linkageRow)
			{
			try
				{
				/*
				System.err.println(getPairOfIndividualIds());
				System.err.println(this.typesAA);
				System.err.println(this.pairid2theoritical);
				System.err.println(this.pairid2couple);
				*/
				if(getPairOfIndividualIds().isEmpty())
					{
					return null;
					}
				
				
				double zScoreAA=0;
				for(Float catIbd0: this.typesAA.keySet())
					{
					int eventsAA=0;
					int totalAA=0;
					for(Integer pair_id:getPairOfIndividualIds())
						{
						if(this.pairid2theoritical.get(pair_id).getIbd0()!=catIbd0) continue;
						
						if(this.pairid2couple.get(pair_id).getConcordance()!=PairOfIndividual.Concordance.AFFECTED)
							{
							continue;
							}
						totalAA++;
						
						if(getIbd0(linkageRow,pair_id)!=1) continue;
						eventsAA++;
						}
					
					double prAA = pbinom( eventsAA-1, totalAA, catIbd0,true,false ) +
							dbinom(eventsAA,totalAA, catIbd0,false)
							; 
					
					prAA=(prAA>=1.0?0.0999999:prAA);
					
					zScoreAA +=
							qnorm(prAA,0,1.0,false,false)*Math.sqrt(totalAA);
					
					}
				double saveZscoreAA=zScoreAA;
				zScoreAA=zScoreAA/Math.sqrt(nbPairsAA);//truc a utiliser LA ICI 
				
				
				
				
				double zScoreAU=0;
				for(Float catIbd0: this.typesAU.keySet())
					{
					int eventsAU=0;
					int totalAU=0;
					for(Integer pair_id: getPairOfIndividualIds())
						{
						if(this.pairid2theoritical.get(pair_id).getIbd0()!=catIbd0) continue;
						
						if(this.pairid2couple.get(pair_id).getConcordance()!=PairOfIndividual.Concordance.DISCORDANT) continue;
						totalAU++;
						
						if(getIbd0(linkageRow,pair_id)!=1) continue;
						eventsAU++;
						}
					
					double prAU = pbinom( eventsAU-1, totalAU, catIbd0,false,false) +
							dbinom(eventsAU,totalAU, catIbd0,false)
							; 
					
					prAU=(prAU>=1.0?0.0999999:prAU);
					
					zScoreAU +=
							qnorm(prAU,0,1.0,false,false)*Math.sqrt(totalAU);
					
					}
				double saveZscoreAU=zScoreAU;
				zScoreAU=zScoreAU/Math.sqrt(nbPairsAU);
	
				Stat stat=new Stat();
				stat.zScoreAA=zScoreAA;
				stat.zScore= (saveZscoreAA+saveZscoreAU)/Math.sqrt(nbPairsAA+nbPairsAU);
				
				return stat;
				 } 
			catch (Exception e) {
				e.printStackTrace();
			    return null;
			    }
			catch (Throwable e) {
				e.printStackTrace();
			    return null;
			    }
			}
		@Override
		public Double getValue(LinkageRow linkageRow)
			{
			
			Stat stat=getStats(linkageRow);
			if(stat==null || stat.isNan())
				{
				return null;
				}
			
			Double v= getValue(stat);
			return Double.isNaN(v) || Double.isInfinite(v) ? null:v;
			}
		protected abstract Double getValue(Stat s);
		}
	
	private class DataExtractorKinshipZScore
		extends DataExtractorKinship
			{
			@Override
			protected Double getValue(Stat stat)
				{
				return stat.lodScore();
				}
			}
	
	private class DataExtractorKinshipZScoreAA
	extends DataExtractorKinship
		{
		@Override
		protected Double getValue(Stat stat)
			{
			return stat.lodScoreAA();
			}
		}
	
	/** Constructor */
	public DataExtractorFactoryIKinship(boolean use_zScoreAA)
		{
		this.use_zScoreAA=use_zScoreAA;
		}
	
	/** Constructor */
	public DataExtractorFactoryIKinship(DataExtractorFactoryIKinship cp)
		{
		super(cp);
		this.use_zScoreAA=cp.use_zScoreAA;
		}
	
	
	@Override
	public Object clone()
		{
		return new DataExtractorFactoryIKinship(this);
		}
	
	
	@Override
	public String getName()
		{
		return "Kinship: "+(this.use_zScoreAA?"Lod Score AA":"Lod Score");
		}
	
	@Override
	public String getDescription() {
		return getName()+" p:="+getTreshold();
		}
	
	public DataExtractor createDataExtractor()
		{
		DataExtractorKinship ex=null;
		if(use_zScoreAA)
			{
			ex=new DataExtractorKinshipZScoreAA();
			}
		else
			{
			ex=new DataExtractorKinshipZScore();
			}
		ex.setModel(this.getModel());
		ex.setTransaction(this.getTransaction());
		ex.setTreshold(this.getTreshold());
		return ex;
		}
	
	
	
	
	
	
	
	
	
	
	
	
	
	

	  /**
	   *
	   * Calculates the value of the cumulative distribution function
	   *
	   * @param dist the distribution
	   * @param q the value
	   * @param lowerTail if true, return the value P(x < q), otherwise P(x > q)
	   * @param logP  if true, return the natural logarithm of the probability
	   * @return  the probability that the random variable will take the value less than (greater than)
	   * {@code q}
	   */
	  private static double _p(Distribution dist, double q, boolean lowerTail, boolean logP) {
	    double p;
	    try {
	      p = dist.cumulativeProbability(q);
	    } catch (MathException e) {
	    	e.printStackTrace();
	      return Double.NaN;
	    } catch (MathRuntimeException e) {
	    	e.printStackTrace();
	      return Double.NaN;
	    }
	    if (!lowerTail) {
	      p = 1.0 - p;
	    }
	    if (logP) {
	      p = Math.log(p);
	    }

	    return p;
	  }

	  private static double _d(IntegerDistribution dist, double x, boolean log) {
		    double d = dist.probability(x);
		    if (log) {
		      d = Math.log(d);
		    }
		    return d;
		  }


	  private static double _q(ContinuousDistribution dist, double p, boolean lowerTail, boolean logP) {
		    if (logP) {
		      p = Math.exp(p);
		    }
		    double q = 0;
		    try {
		      q = dist.inverseCumulativeProbability(p);
		    } catch (IllegalArgumentException e) {
		    	e.printStackTrace();
		      return Double.NaN;
		    } catch (MathException e) {
		    	e.printStackTrace();
		      return Double.NaN;
		    } catch (MathRuntimeException e) {
		    	e.printStackTrace();
		      return Double.NaN;
		    }
		    if (!lowerTail) {
		      q = -q;
		    }
		    return q;
		  }


	
	private static double pbinom (int x, int size, double prob,boolean lowerTail,boolean logP) 
		{
		 return _p(new BinomialDistributionImpl(size, prob), x, lowerTail, logP);
		}
	
	private static double  dbinom(double x,  int size, double prob, boolean log)
		{
		return _d(new BinomialDistributionImpl(size, prob), x, log);
		}
		
	private static double qnorm(double p,double mean, double sd, boolean lowerTail, boolean logP)
		{
		return _q(new NormalDistributionImpl(mean, sd), p, lowerTail, logP);
		}
	}
